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Abstract. The Wilson fermion determinant can be written as product of the de- 
terminants of two hermitian positive definite matrices. This formulation allows to 
simulate non-degenerate quark flavors by means of the hybrid Monte Carlo al- 
gorithm. A major numerical difficulty is the occurrence of nested inversions. We 
construct a Uzawa iteration scheme which treats the nested system within one 
iterative process. 



1 Introduction 

The strong interaction between the quarks is described by quantum chromo- 
dynamics (QCD), a constitutive part of the standard model of elementary 
particle physics. QCD develops a strong coupling constant at the hadronic 
mass scale and thus, perturbation theory cannot be applied. The only known 
non-perturbative ab-initio method is to simulate QCD on a 4-dimensional 
space-time lattice by use of Monte Carlo methods as known from statistical 
physics. 

One would think that lattice simulations of QCD have to take into ac- 
count virtual loops from all six quark flavors, up (u), down (d), strange (s), 
charm (c) bottom (b), and top (t). Their masses span a wide range from 
about 3 MeV to about 180 GeV. 1 The masses of the three heavy quarks 
lie above the momentum scale set by dynamical chiral symmetry breaking. 
The QCD Lagrangian is chirally symmetric for vanishing quark masses. Chi- 
ral symmetry is explicitly broken by the quark masses and spontaneously 
by the dynamics. For light quarks, (u, d, s), dynamical breaking dominates, 
as signaled by zero mass Goldstone bosons, i.e., the pion octet. Supposedly 
only the light quarks contribute with virtual loops; the contributions of the 
heavier ones are assumed to be negligible. 

1 The light quark (u, d, s) masses are so-called "current" masses determined within 
the MS scheme at a renormalization scale = 2 GeV. The masses of the heavy 
quarks (c, b, t) are "running" masses determined at /i = m c ^ t in the MS scheme 

[!]• 
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A straightforward lattice QCD simulation including u, d, and s is difficult: 
on one hand, u and d are very light. As we know from chiral perturbation 
theory, the square of the pion mass is proportional to the quark mass [2]. 
Therefore, the pion acquires a small mass, i.e. a large correlation length, 
£ = mo' f° r u an ^ ^ masses approaching the chiral limit, m u and = 0. a 
denotes the lattice spacing. The lattice volume must be larger by more than 
one order of magnitude than possible today to reach physical parameter val- 
ues for m u and m^. Hence, one has to recourse to extrapolations using results 
of simulations at artificial parameters for m u and and small lattices, far 
off the chiral limit. These difficulties go along with an increasing condition 
number of the fermionic matrix M, therefore its inversion suffers from critical 
slowing down approaching the chiral limit [3,4]. 

On the other hand, until recently, the only exact simulation algorithm 
for QCD with dynamical fermions 2 was the hybrid Monte Carlo algorithm 
(HMC) [5,6]. The benefits of HMC in reducing the computational complexity 
are achieved by treating the fermionic determinants as a stochastic estimate 
using Gaussian random fields [7] , with the advantage, that instead of comput- 
ing the determinant, only the solution of a linear system is required. For this 
approach the fermion matrix M must be hermitian positive definite (hpd). 

The Wilson fermion discretization of the fermionic sector of QCD de- 
scribes single flavors [8] , in contrast to staggered fermions that represent four 
fermions, intermixed in spin-flavor space. 3 Unfortunately, Wilson fermions 
must violate chiral symmetry, see Ref. [9]. As a consequence, the Wilson 
fermion matrix is complex non-normal and thus cannot be included in the 
HMC scheme. To avoid this problem, one usually simulates two mass degen- 
erate light quarks, an approximation justified by the fact that both quarks, 
u and d, are light and close in mass compared to the next heavier one, s. 
The product of two mass degenerate determinants indeed amounts to an hpd 
matrix, that can be simulated by the HMC. The minor price to pay is that 
the two light quarks are mass-degenerate, the major price is that the s quark 
is not included in the simulation. 

There are attempts to evaluate operators, that contain the s quark, by 
application of a "partially quenched analysis" [10]. The PQA extrapolates 
only one of the valence quarks towards the chiral limit and holds the other 
one at the mass of the s quark — of course within the light mass-degenerate 
u-d-sea. PQA, however, leads to inconsistencies: the so-called J-parameter 
does not acquire its physical value [11,12]. 

One can try to employ an approximate one-flavor determinant represen- 
tation. Such simulations, however, are plagued by systematic uncertainties, 
as one has to recourse to non-exact simulation algorithms like the hybrid 
molecular dynamics algorithm [13,14]. 

2 "Dynamical" in contrast to "quenched" simulations that neglect fermionic loops. 

3 The naive discretization of the Dirac operator yields, in addition to the true 
mode, 15 "doublers". 
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On the other hand, in the framework of the multi-boson algorithm [15], 
one can define an exact single flavor simulation scheme for Wilson fermions 
[16-18]. Here many additional bosonic auxiliary fields are required, and it is 
not clear if the method would be well suited for more complicated actions or 
allows to exploit sophisticated preconditioning. 

In this paper, we propose to include the Wilson fermion determinant into 
HMC such that a single quark can be described. We show that the fermionic 
determinant can be written as a product of the determinants of two hpd 
matrices. This is achieved by the representation of the fermion matrix via its 
Schur complement. Both determinants can be treated by the hybrid Monte 
Carlo in the standard manner. There are, however, some caveats: one of 
the two matrices involves the inverse of the other one. As a consequence, a 
nested inversion must be carried out. We propose to use a Uzawa-like inverter 
which allows to solve the nested system within one iterative process. A second 
difficulty stems from the initialization of the fermion action in the HMC 
algorithm: since we have a single flavor system, we have to compute a square 
root at the beginning of a HMC trajectory. 

The paper is organized as follows: in section 2 we introduce the basics of 
lattice QCD and will discuss the difficulties simulating single flavors of Wilson 
fermions by the HMC. In section 3 we give the transformation of the fermion 
matrix using the Schur complement, in section 4 we present the non-mass- 
degenerate HMC for Wilson fermions along with a discussion of numerical 
problems. In section 5, we introduce an algorithm that can compute the 
nested inversion in one iterative process. Finally, we summarize and give an 
outlook. 

2 Numerical Problem 

The action of lattice QCD is the sum of a gauge part (which is not relevant 
for the following) and a fermionic part: 

S[U,i>,^}:=S g [U}+ J2 (1) 

x,y,a,b,a,/3 

V is 4-dimensional volume in Euclidean space. The full matrix is a tensor 
product with 3x3 SU(3) matrices U (color) and four 4x4 matrices 7 M 
(Dirac spin), hence 

MeC 12Vxl2V (2) 

Here the indices x an y stands for the 4-dimcnsional coordinates, a and b 
denote color space and a and (3 Dirac space indices. 

The fermion fields ip x are Grassmann variables. Therefore, the Grass- 
mann integral of the exponential of the fermionic part of the action, i.e., the 
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fermionic part of the path integral, can be carried out, 

[#][#] e - * M * =det(M), (3) 



/> 



to yield the determinant of M for arbitrary complex M. 
Wilson fermions are defined by the interaction matrix 

M = 1-kD, MeC 12Vxl2V , 

4 

Ac,y = (1 - 7m) M x ) ^,y-„ + (1 + 7m) (X - AO ftc,y+M- ( 4 ) 

The Dirac matrices in Euclidean space satisfy the anti-commutator rela- 
tions 

7(u7^ + 7^7(u = 2 <W ( 5 ) 

It is important for the following considerations to work with the chiral rep- 
resentation of the Dirac matrices, as given in the appendix A. In this repre- 
sentation, 75, defined as the product 75 = 74717273, is diagonal. 

M is a complex non-hermitian matrix which moreover is not normal, i.e., 

MM* ^ M*M. (6) 

Thus M cannot be diagonalized by a unitary transformation. At most, it 
might be diagonalizable by a similarity transformation. However, M exhibits 
the so-called 75-symmetry or 75-hermiticity: 

M75 = 75M. (7) 

Therefore, the matrix 

Q ■■= 75M (8) 

is hermitian. We can immediately read off the chiral representation of the 
Dirac matrices, see the appendix, that, for k = 0, Q exhibits an equal number 
of positive and negative eigenvalues. In general, Q is indefinite. 

The matrix M^M is hpd. Since the determinant of M^M represents two 
mass degenerate fermion flavors, 

det (M t Af) = det (Af ) det (Af) , (9) 

the two flavor situation is ideal for HMC simulations. Being unitary diag- 
onalizable it can be represented by a Gaussian integral over complex fields 



det (Af) det (M) = (^j^ J [#] [d(f>*]e~^ ( MtM ) (i ) 
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For the single flavor M, such a simple construction fails even for the 
situation n < k c where M is positive definite, because the Jacobian of the 
similarity transform — if the diagonalization is feasible at all — is not equal to 
1. 

On the other hand, 

det(M) =det(Q). (11) 

But, Q is indefinite. Therefore, a Gaussian integral representation for its 
determinant does not exist. 

We conclude that so far there is no direct way to include individual flavors 
within the HMC algorithm. 



3 Schur Complement of Q 



Let's consider the hermitian matrix Q with the Dirac matrices given in the 
chiral representation of appendix A. The explicit form of Q in the chiral 
representation is given by 



Q 



Qn Q12 
Q21 Q22 



I(I-kDh) -kD 12 
t 

12 



-kD{ 2 -1(1-kD u )J' 



(12) 



Here, the bold face 1 represents a unit matrix in 2 x 2 spin space. D\\ does 
not carry spin indices while D 12 consists of 2 x 2 blocks in spin space, 



4 

D12 = ^2v» [PpM ^x,x+ M - U£(x - /z) 5 x ,x- M ] , 



with 



Vl = icr li V2 = «CT 2 , 



(13) 



(14) 



The Schur complement of a 2 x 2 block matrix with non-singular block A 
is given by 



'A B~ 




1 0" 




A 




1 A- X B 


CD 




CA' 1 1 




D-CA^B 




1 



(15) 



Hence 



det(Q) = det (l(l - kD u )) det (l(l - kD u ) + k 2 D 12 [1{1 - kAi)]" 1 !'^ • 

(16) 
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The minus sign of the Q22 term cancels as the rank of the matrix is even. From 
Eq. 13 we know that the first matrix is diagonal in Dirac space. Therefore its 
determinant can be written as follows: 

det(Q) = dct (1(1 - k£>h)) = [det(l - kDu)} 2 = det(l - k£>ii) 2 , (17) 

where Q\ = (1 — Dn) 2 is hpd for any n. Q w carries no spin index. 
The second matrix 

Q sc = 1(1 - kD u ) + k 2 D 12 [1(1 - kD 11 )]- 1 d\ 2i (18) 

the Schur complement, is hermitian, as it is the sum of hermitian matrices. 
For k = 0, it is hpd, thus there exists some value of k, k^ c , for which it 
becomes indefinite. For k < k^ c , Q sc is hpd. 

4 One Flavor HMC 

We have shown that 

dct (M) = dct (Q) = det (Q w ) dct (Q sc ), (19) 

for Qu being non-singular. If k < k sc , both matrices are hpd and can be 
represented by Gaussian integrals. Let (p € C 3V and \ G C 6V : 

dct (Q w ) dct (Q sc ) 

= (^) 3y (^) 6V |[#][#*][dx][dx1e-* t W-)- a ^ t W«)- 1 x. 

(20) 

The fermionic part of the path integral being expressed in terms of the 
pseudo-fcrmion degrees of freedom 4> and %, the HMC algorithm can proceed 
as usual [19], see the ^-algorithm in Ref. [6]. The heat bath for the (f> fields 
is trivial like in the previous case of mass-degenerate QCD. Let R e C 3V be 
Gaussian distributed. With 

4 = QtR, (21) 

R}R = 4>\Q w )~ 2 (t>. (22) 

and 

{2nf v det{Ql) = J [d4>]W]e~^ . (23) 

The refreshment of the \ fields causes more trouble, however. In order to 
generate the required distribution, we have to solve a linear system involving 
the square root of Q sc . Let R € C 6V be Gaussian distributed: 



x = Q sc l2 R- 



(24) 
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Then 

rtR = xHQsc)- 1 X (25) 

and 

(27r) 6V dct(Q sc )= J[d X ][dx*]e-^ Q -J~ lx (26) 

Even though the solution of Eq. 25 is only required at the beginning of a 
trajectory, it is particularly expensive as we must cope with the internal 
inversion of Q\\. The implications of this issue arc not clarified yet. As an 
alternative, we might include the dynamics of the X fields into HMC avoiding 
the expensive heat bath step for Q sc , similar to the i?-algorithm of Ref. [6]. 



5 Inexact Uzawa Iterations 

The nested inversion within the solution of 

QscX = X (27) 

has to be carried out at each molecular dynamics time step of HMC as well 
as in the computation of the action. By Uzawa iterations [20] we avoid the 
solution of the linear system represented by Q w within each iteration step of 
the Q sc iteration. 

The problem to solve is given by: 

(Q w + D 12 Q- 1 D\ 2 )X = X , (28) 

with x,i>e C 6V . 

Let's consider the Jacobi iteration for Q sc , 

* (<+1) = X + nD n X^ - Dia Q- 1 D\ 2 X®. (29) 

The fixed-point solution X of Eq. 29, defined by 

||j5f( i+1 ) -X®\\ < e, e^O (30) 

solves equation Eq. 28, as can be easily verified. 

Next we transform the simple Jacobi iteration into the Uzawa iteration: 

X {l+1) = x + kDuX® - £>ia r w 

Y ii+1) =Qw lD i2 x{i) - (31) 

At this stage, the Uzawa scheme would certainly not offer an advantage as still 
we need to carry out an inversion in each iteration step. Let's turn towards 
the inexact Uzawa iteration. Choose 



(32) 
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as a preconditioner for Q w which is easily invertible. 

X il+1) = X + kDuX® - D w Y {1) 

= yW + p-i(Dt 2 A-W - Y®) (33) 

The fixed-point solution of Eq. 33, defined by 

-x®\\ < e , 
||y(i+i) _y(0|| < e '. (34) 

for e, e' — > 0, solves Eq. 28. Using the second relation, one finds 

o = p- 1 (-Q tu r + J Dt 2 x), 

oy = Q- 1 i?l 2 X (35) 
As a cheap preconditioner we can use a truncated Neumann series for Q" 1 : 

P- 1 = (1- KD^y 1 = 1 + kD u + ... (36) 

One has to find the optimal length of this series compared to the convergence 
of the outer iteration. In our implementation we have chosen the simple first- 
order approximation. 

P- 1 = 1 + kDu. (37) 

We have made first tests of the Uzawa iteration in comparison to a con- 
jugate gradient (CG) for the inversion of Eq. 28. We used a quenched 16 4 
configuration at = 6.0 at a moderate K-value of k = 0.12. The outer CG 
took about 240 iterations while the inner iteration took about 50 steps. The 
Uzawa iteration required about 610 steps in comparison. The overall improve- 
ment for this setting is about a factor of 10. Of course, these results are very 
preliminary and have to be confirmed for other values of k and with real- 
istic HMC simulations. It would also be desirable to implement the Uzawa 
scheme within Krylov subspace solvers. 

For free theory (U = 1), one can diagonalize Dn by Fourier transfor- 
mation. As in the case of Wilson fermions, one finds k c = |. We have, so 
far, investigated the lowest lying eigenvalues of Q w for the interacting case: 
on a 16 3 x 32 configuration, taken from the SESAM ensembles at = 5.6 
and Ksea = 0.1575, we found for the inverse of the critical eigenvalue of 
Du, 1/Xmin — 0.144. The critical k of Dn is reached for a K-value being 
smaller than the critical k value of the Wilson matrix. Therefore, the singu- 
larity might be a barrier screening the approach to the chiral limit. These 
questions have to be clarified in an actual HMC simulation with one-flavor 
QCD. 
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6 Summary and Outlook 

The determinant of the Wilson fermion matrix is equivalent to the product of 
the determinants of two hermitian positive definite matrices. This formulation 
allows to simulate quark flavors with individual mass values by means of the 
HMC algorithm. The solution of the ensuing linear system is hampered by 
a nested inversion problem. However, it can be treated by use of Uzawa 
iterations. 

Currently we perform a feasibility study to gain first experiences with the 
new method. In particular, we are interested in the analytic structure of Q w 
and Q sc . Furthermore, we construct an efficient heat-bath for the refreshment 
of the pseudo-fermions for Q sc based on the polynomial representation of the 
square root of a truncated polynomial that represents Q sc . 

We have presented the method for standard Wilson fermions. The Schur 
decomposition can as well be applied to improved Wilson fermions with clover 
term. Work on this line is in progress. 
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A Chiral Dirac Matrices 



Euclidean Dirac matrices in the chiral representation: 



7i = 



72 



73 = 



74 



75 



/0 -i\ 
-iO 
i 
\i 00 
(0 -A 
10 
10 
\-1000 ) 
(00 -i0\ 

i 

1 
\0-i0 0} 

(0 1 0\ 
1 
10 
\0 1 0) 
(\00 
10 
0-10 
\000 -1/ 



-o-i 

(71 



-<7 2 

cr 2 



-<7 3 
(73 



1 

1 



(38) 
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